load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

  cs = "EXP01"
  myfont = 0.016

  path = "../data/"

  load "loadit"

  xx = "_"

  varname = "Q"

  plotform = "pdf"
  plotname = "figure_timeheight_" + varname + xx + cs 

  ;;.......................
  ;; read data
  ;;.......................

  fna = path +cs+"_dt1200.nc"
  fnb = path +cs+"_dt600.nc"
  fnc = path +cs+"_dt300.nc"
  fnd = path +cs+"_dt150.nc"
  fne = path +cs+"_dt75.nc"

  system("ncdump -c "+fna+" > info ") 

  fla = addfile(fna,"r")
  flb = addfile(fnb,"r")
  flc = addfile(fnc,"r")
  fld = addfile(fnd,"r")
  fle = addfile(fne,"r")

  xa = fla->$varname$ ;;NITEND_MTB
  xb = flb->$varname$ ;;NITEND_MTB
  xc = flc->$varname$ ;;NITEND_MTB
  xd = fld->$varname$ ;;NITEND_MTB
  xe = fle->$varname$ ;;NITEND_MTB

  nt   = dimsizes(xa(:,0,0,0))
  nlev = dimsizes(xb(0,:,0,0))

  v1 = xa(:,:,0,0)
  v2 = xb(:,:,0,0)
  v3 = xc(:,:,0,0)
  v4 = xd(:,:,0,0)
  v5 = xe(:,:,0,0)

  vp = new((/5,nlev,nt/),"float")


 wks = gsn_open_wks(plotform,plotname) 
 gsn_define_colormap(wks,"BlueYellowRed") ; choose colormap

; Individual plot resources
  
  res = True
  res@gsnDraw                = False
  res@gsnFrame               = False

  res@cnFillMode = "RasterFill"
 
  res@lbBoxMinorExtentF     = 0.15    
  res@lbOrientation         = "Vertical"
  res@lbLabelsOn            = True
  res@lbLabelFontHeightF    = 0.02
  ;;res@lbLabelAutoStride     = True
  res@lbLabelStride  = 1
  res@lbTitleFontHeightF    = 0.015
  res@tiXAxisFontHeightF    = 0.02
  res@tiXAxisString         = "Days" 
  res@tiYAxisFontHeightF    = 0.02
  res@tiYAxisString         = "Pressure (mb)"
   
;  res@tmXBMinorOn           = False
  res@tmYLMinorOn           = False
  res@tmXTMinorOn           = False
  res@tmYRMinorOn           = False
  res@tmXTOn                = False
  
  ;;res@gsnRightString        = "" 
  ;;res@gsnLeftString         = "" 
  res@gsnStringFontHeightF  = 0.02
  res@gsnCenterStringFontHeightF  = 0.03
  
  res@cnFillOn             = True          ; turn on color
  res@gsnSpreadColors      = True     ; use full colormap
  res@cnLinesOn            = False    ; no contour lines
 
  res@gsnSpreadColors      = True          ; use full range of colormap
  
  ;;if (cmin.ne.cmax) then 
  ;;  res@cnLevelSelectionMode = "ManualLevels"
  ;;  res@cnLevelSpacingF      = cgap
  ;;  res@cnMinLevelValF       = cmin
  ;;  res@cnMaxLevelValF       = cmax
  ;;end if
  res@vpWidthF             = 0.8          ; change aspect ratio of plot
  res@vpHeightF            = 0.3
  ;;res@trYReverse           = True
  
  v2!0 = "time"  ;Reconstruct meta data as the time dimension has changed.
  v2!1 = "lev"
  
  v2&time = xa&time
  v2&lev  = xa&lev
  v2&lev@units = "hPa"
  v2&lev@long_name = "Pressure Levels"

  
  v1!0 = "time"
  v1!1 = "lev"
  
  v1&time = xa&time
  v1&lev  = xa&lev
  v1&lev@units = "hPa"
  v1&lev@long_name = "Pressure Levels"

  v3!0 = "time"
  v3!1 = "lev"

  v3&time = xa&time
  v3&lev  = xa&lev
  v3&lev@units = "hPa"
  v3&lev@long_name = "Pressure Levels"

  v4!0 = "time"
  v4!1 = "lev"

  v4&time = xa&time
  v4&lev  = xa&lev
  v4&lev@units = "hPa"
  v4&lev@long_name = "Pressure Levels"

  v5!0 = "time"
  v5!1 = "lev"

  v5&time = xa&time
  v5&lev  = xa&lev
  v5&lev@units = "hPa"
  v5&lev@long_name = "Pressure Levels"


  res@tmYRMode             = "Automatic"          ; turn off special labels on right axis

  plot = new(5,graphic)   

  res@cnLevelSelectionMode = "ExplicitLevels"    			; set explicit contour levels
  res@cnLevels    = (/ 1.e-5,2.e-5,5.e-5,1.e-4,2.e-4,5.e-4,1.e-3,2.e-3,5.e-3,1.e-2/) 

  res@gsnLeftString = "a) dt1200" 
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~" 
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~" 

  plot(0) = gsn_csm_pres_hgt(wks,v1({lev|100:1000},time|:),res)

  ;;v1 = where(v1.lt.100.,-999,v1) 
  ;;v2 = where(v2.lt.100.,-999,v2) 
  ;;v3 = where(v3.lt.100.,-999,v3) 
  ;;v4 = where(v4.lt.100.,-999,v4) 
  ;;v5 = where(v5.lt.100.,-999,v5) 

  v1@_FillValue = -999  
  v2@_FillValue = -999  
  v3@_FillValue = -999  
  v4@_FillValue = -999  
  v5@_FillValue = -999  

  v2 = v2 - v1 
  v3 = v3 - v1 
  v4 = v4 - v1 
  v5 = v5 - v1 

  delete(res@cnLevels) 
  ;;delete(res@cnLevelSelectionMode) 
  res@cnLevels    = (/ -2.e-3,-1.e-3,-5.e-4,-2.e-4,-1.e-4,-5.e-5,-2.e-5,-1.e-5,0.,1.e-5,2.e-5,5.e-5,1.e-4,2.e-4,5.e-4,1.e-3,2.e-3/) 

  res@gsnLeftString = "b) dt600"  
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~" 
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~" 
  plot(1) = gsn_csm_pres_hgt(wks,v2({lev|100:1000},time|:),res)

  res@gsnLeftString = "c) dt300" 
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~" 
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~" 
  plot(2) = gsn_csm_pres_hgt(wks,v3({lev|100:1000},time|:),res)

  res@gsnLeftString = "d) dt150"
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~"
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~"
  plot(3) = gsn_csm_pres_hgt(wks,v4({lev|100:1000},time|:),res)

  res@gsnLeftString = "e) dt75"
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~"
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~"
  plot(4) = gsn_csm_pres_hgt(wks,v5({lev|100:1000},time|:),res)


;; draw panel with white space added
 resP                 = True
 resP@gsnPanelYWhiteSpacePercent = 5
 resP@gsnPanelXWhiteSpacePercent = 5
 gsn_panel(wks,plot,(/5,1/),resP)


end







